{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 结构能量打分\n",
    "\n",
    "## 概述\n",
    "\n",
    "rosetta使用基于物理和统计的能量函数项,计算生物大分子的能量.本章结束应当理解:\n",
    "\n",
    "- rosetta中计算生物大分子的能量意味着什么\n",
    "- 一个特定的相互作用对总能量/分数的贡献\n",
    "- 如何解释和比较能量分数\n",
    "- 如何准备用于打分的生物大分子的结构\n",
    "- 如何对生物大分子打分\n",
    "- 如何更用其他能量函数\n",
    "- 如何自定义打分函数的能量项\n",
    "- 如何获得能量分数的每个残基贡献\n",
    "\n",
    "## rosetta的打分\n",
    "\n",
    "在rosetta中,打分scoring即是在计算生物大分子的能量.\n",
    "\n",
    "rosetta有一个打分函数ref2015,用于计算L型氨基酸组成的蛋白的全原子相互作用能量.也存在其他特定应用的全原子打分函数和质心打分函数.你还可以自定义打分函数.\n",
    "\n",
    "## 打分函数\n",
    "\n",
    "打分函数是能量项的加权求和,有些代表物理相互作用(静电力,vdW等),其他代表统计项(如Ramachandran空间中的扭转角概率).以下是ref2015中的能量项.\n",
    "|能量项|描述|\n",
    "|-|-|\n",
    "|fa_atr|不同残基原子间的Lennard-Jones吸引作用|\n",
    "|fa_rep|不同残基原子间的Lennard-Jones排斥作用|\n",
    "|fa_sol|Lazaridis-Karplus溶剂化能|\n",
    "|fa_intra_sol_xover4|残基内Lazaridis-Karplus溶剂化能|\n",
    "|lk_ball_wtd|不对称溶剂化能|\n",
    "|fa_intra_rep|相同残基原子间Lennard-Jones排斥作用|\n",
    "|fa_elec|与距离相关的库仑静电势|\n",
    "|pro_close|脯氨酸闭环能量,以及前驱残基的psi角能量|\n",
    "|hbond_sr_bb|一级序列中的骨架-骨架近端氢键|\n",
    "|hbond_lr_bb|一级序列中的骨架-骨架远端氢键|\n",
    "|hbond_bb_sc|侧链-骨架氢键|\n",
    "|hbond_sc|侧链-侧链氢键|\n",
    "|dslf_fa13|双硫键几何势|\n",
    "|rama_prepro|Ramachandran倾向(脯氨酸和其他残基的查询表格分开)|\n",
    "|omega|骨架中的Omega二面角. 是平面上的谐波约束,标准差约为6度|\n",
    "|p_aa_pp|给定phi和psi扭转角的氨基酸概率|\n",
    "|fa_dun|由Dunbrack统计库衍生的侧链旋转异构内能|\n",
    "|yhh_planarity|保证酪氨酸羟基在芳香环平面上的特殊扭转角势能|\n",
    "|ref|每个氨基酸的参考能量.平衡了氨基酸内能项,主要在设计中起作用|\n",
    "|METHOD_WEIGHTS|本身不是能量,而是每个氨基酸的参数,被ref能量项使用|\n",
    "\n",
    "更进一步的描述参见[分数类型](#TODO).\n",
    "\n",
    "能量项权重如下:\n",
    "\n",
    "```txt\n",
    "METHOD_WEIGHTS ref 0.773742 0.443793 -1.63002 -1.96094 0.61937 0.173326 0.388298 1.0806 -0.358574 0.761128 0.249477 -1.19118 -0.250485 -1.51717 -0.32436 0.165383 0.20134 0.979644 1.23413 0.162496 \n",
    "fa_atr 1\n",
    "fa_rep 0.55\n",
    "fa_sol 0.9375\n",
    "fa_intra_rep 0.005\n",
    "fa_elec 0.875\n",
    "pro_close 1.25\n",
    "hbond_sr_bb 1.17\n",
    "hbond_lr_bb 1.17\n",
    "hbond_bb_sc 1.17\n",
    "hbond_sc 1.1\n",
    "dslf_fa13 1.25\n",
    "rama 0.25\n",
    "omega 0.625\n",
    "fa_dun 0.7\n",
    "p_aa_pp 0.4\n",
    "yhh_planarity 0.625\n",
    "ref 1\n",
    "```\n",
    "\n",
    "打分算法参见[打分解释](#TODO).\n",
    "\n",
    "## rosetta分数与真实能量的比较\n",
    "\n",
    "尽管许多能量项是基于物理的,但有些统计项是基于已知结构的观测.\n",
    "\n",
    "打分分数越低,更可能接近原生结构,但分数并不能直接转化为物理能量单位(如kcal/mol).我们使用Rosetta单位能量(Rosetta Energy Units,REU)进行表示.\n",
    "\n",
    "译注:对于ref2015打分函数,现在REU已经矫正,近似kcal/mol.\n",
    "\n",
    "同时,由于分数与打分函数有关,通常使用不同的打分函数比较结构的分数是没有意义的.\n",
    "\n",
    "## 示例"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "sys.path.append(os.path.abspath('..'))\n",
    "from note_utils.path import chdir, add_rosetta_env\n",
    "\n",
    "rosetta_home=add_rosetta_env('~/git_develop/rosetta') # rosetta主目录路径自行修改\n",
    "tutorial_path=rosetta_home/'demos/tutorials/scoring'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 准备要打分的PDB文件\n",
    "\n",
    "使用`score_jd2`程序打分,接受PDB输入格式.直接从PDB数据库下载的PDB文件可能不适合直接进行打分,因为输入可能存在各种问题.\n",
    "\n",
    "例如尝试对3TDM结构进行打分:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cd /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "当前目录从 /home/regen/git_develop/rosetta_notebook/00intro \n",
      "临时切换到 /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "core.chemical.GlobalResidueTypeSet: Loading (but possibly not actually using) 'PO4' from the PDB components dictionary for residue type 'pdb_PO4'\n",
      "score.sc\n",
      "当前目录切换回到 /home/regen/git_develop/rosetta_notebook/00intro\n"
     ]
    }
   ],
   "source": [
    "with chdir(tutorial_path):\n",
    "    !$rosetta_home/source/bin/score_jd2.linuxgccrelease -in:file:s input_files/from_rcsb/3tdm.pdb |grep PO4\n",
    "    !ls score.sc\n",
    "    !rm score.sc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于未知的磷酸基团,无法确保对其打分(旧版将直接报错退出).\n",
    "\n",
    "要对该PDB打分,添加`-ignore_unrecognized_res`选项将无视PDB中不认识的基团."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cd /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "当前目录从 /home/regen/git_develop/rosetta_notebook/00intro \n",
      "临时切换到 /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "core.chemical.GlobalResidueTypeSet: Loading (but possibly not actually using) 'PO4' from the PDB components dictionary for residue type 'pdb_PO4'\n",
      "score.sc\n",
      "当前目录切换回到 /home/regen/git_develop/rosetta_notebook/00intro\n"
     ]
    }
   ],
   "source": [
    "with chdir(tutorial_path):\n",
    "    !$rosetta_home/source/bin/score_jd2.linuxgccrelease -in:file:s input_files/from_rcsb/3tdm.pdb -ignore_unrecognized_res | grep PO4\n",
    "    !ls score.sc\n",
    "    !rm score.sc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "该选项在旧版将保证打分应用的正常运行,而在目前版本与默认行为一致.\n",
    "\n",
    "该命令将生成`score.sc`分数文件,要进行下一步,执行`rm score.sc`删除分数文件,否则后续执行的所有打分命令的结果都会附加到该文件中.\n",
    "\n",
    "如果输入PDB中缺失重原子将自动补全.可以指定`-out:pdb`将处理后的结构进行输出.以1QYS为例:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cd /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "当前目录从 /home/regen/git_develop/rosetta_notebook/00intro \n",
      "临时切换到 /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "core.io.pose_from_sfr.PoseFromSFRBuilder: Reading MSE as MET!\n",
      "core.io.pose_from_sfr.PoseFromSFRBuilder: Reading MSE as MET!\n",
      "core.io.pose_from_sfr.PoseFromSFRBuilder: Reading MSE as MET!\n",
      "core.io.pose_from_sfr.PoseFromSFRBuilder: Reading MSE as MET!\n",
      "core.io.pose_from_sfr.PoseFromSFRBuilder: Reading MSE as MET!\n",
      "core.io.pose_from_sfr.PoseFromSFRBuilder: Reading MSE as MET!\n",
      "core.io.pose_from_sfr.PoseFromSFRBuilder: Reading MSE as MET!\n",
      "core.io.pose_from_sfr.PoseFromSFRBuilder: Reading Selenium SE from MSE as SD from MET\n",
      "core.io.pose_from_sfr.PoseFromSFRBuilder: Reading MSE as MET!\n",
      "core.pack.pack_missing_sidechains: packing residue number 13 because of missing atom number 6 atom name  CG \n",
      "core.pack.pack_missing_sidechains: packing residue number 15 because of missing atom number 6 atom name  CG \n",
      "core.pack.pack_missing_sidechains: packing residue number 25 because of missing atom number 6 atom name  OG \n",
      "core.pack.pack_missing_sidechains: packing residue number 28 because of missing atom number 6 atom name  CG \n",
      "core.pack.pack_missing_sidechains: packing residue number 29 because of missing atom number 6 atom name  CG \n",
      "core.pack.pack_missing_sidechains: packing residue number 32 because of missing atom number 6 atom name  CG \n",
      "core.pack.pack_missing_sidechains: packing residue number 34 because of missing atom number 6 atom name  CG \n",
      "core.pack.pack_missing_sidechains: packing residue number 44 because of missing atom number 6 atom name  CG \n",
      "core.pack.pack_missing_sidechains: packing residue number 45 because of missing atom number 6 atom name  CG \n",
      "core.pack.pack_missing_sidechains: packing residue number 53 because of missing atom number 6 atom name  CG \n",
      "core.pack.pack_missing_sidechains: packing residue number 60 because of missing atom number 6 atom name  CG \n",
      "core.pack.pack_missing_sidechains: packing residue number 71 because of missing atom number 6 atom name  CG \n",
      "1qys_0001.pdb  score.sc\n",
      "当前目录切换回到 /home/regen/git_develop/rosetta_notebook/00intro\n"
     ]
    }
   ],
   "source": [
    "with chdir(tutorial_path):\n",
    "    !$rosetta_home/source/bin/score_jd2.linuxgccrelease -in:file:s input_files/from_rcsb/1qys.pdb -out:pdb -overwrite | grep -E 'MSE|packing'\n",
    "    !ls score.sc 1qys*.pdb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "rosetta会将MSE转换为MET,补全缺失原子.将生成`score.sc`和`1qys_0001.pdb`文件.注意补全操作是随机的,每次补全的结构和分数文件默认有微小区别.\n",
    "\n",
    "新生成的PDB文件会包含缺失的原子,并将MSE替换为MET.这个结构可以准确打分.\n",
    "\n",
    "打分文件中`total_score`是很大的正数,说明该结构并不被看好.\n",
    "\n",
    "**这不代表结构不稳定,而是rosetta认为该PDB文件中很可能存在立体冲突**.\n",
    "\n",
    "分数文件的格式类似:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cd /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "当前目录从 /home/regen/git_develop/rosetta_notebook/00intro \n",
      "临时切换到 /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "SEQUENCE: \n",
      "SCORE: total_score dslf_fa13    fa_atr    fa_dun   fa_elec fa_intra_rep fa_intra_sol_xover4              fa_rep              fa_sol hbond_bb_sc hbond_lr_bb    hbond_sc hbond_sr_bb linear_chainbreak lk_ball_wtd       omega overlap_chainbreak            p_aa_pp pro_close rama_prepro         ref yhh_planarity description \n",
      "SCORE:     206.215     0.000  -499.471   298.401  -110.945        1.440              19.595             207.419             316.610      -2.475     -31.013      -3.900     -29.587             0.000      -8.147       0.951              0.000             -8.414     0.000      31.233      24.518         0.000 1qys_0001\n",
      "当前目录切换回到 /home/regen/git_develop/rosetta_notebook/00intro\n"
     ]
    }
   ],
   "source": [
    "with chdir(tutorial_path):\n",
    "    !cat score.sc\n",
    "    !rm score.sc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "推荐总是对PDB先使用驰豫流程再打分,将释放立体冲突.更多驰豫细节见后续教程.\n",
    "\n",
    "### 打分基础\n",
    "\n",
    "本节使用经过优化的1QYS结构进行打分.首先使用默认打分函数ref2015.\n",
    "\n",
    "本节开始我们使用flag文件,不使用空格分隔参数的方式传入选项,因为我们要设置的选项会逐渐变多.\n",
    "\n",
    "目前的flag文件只包含输入文件和输出文件路径:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cd /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "当前目录从 /home/regen/git_develop/rosetta_notebook/00intro \n",
      "临时切换到 /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "-in:file:s input_files/1qys.pdb\n",
      "\n",
      "-out:file:scorefile output_files/score.sc\n",
      "当前目录切换回到 /home/regen/git_develop/rosetta_notebook/00intro\n"
     ]
    }
   ],
   "source": [
    "with chdir(tutorial_path):\n",
    "    !cat flag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用参数文件进行打分的命令如下:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cd /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "当前目录从 /home/regen/git_develop/rosetta_notebook/00intro \n",
      "临时切换到 /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "SEQUENCE: \n",
      "SCORE: total_score       score dslf_fa13    fa_atr    fa_dun   fa_elec fa_intra_rep fa_intra_sol_xover4              fa_rep              fa_sol hbond_bb_sc hbond_lr_bb    hbond_sc hbond_sr_bb linear_chainbreak lk_ball_wtd       omega overlap_chainbreak            p_aa_pp pro_close rama_prepro         ref        time yhh_planarity description \n",
      "SCORE:    -224.872    -224.872     0.000  -501.794   104.484  -146.069        1.118              18.786              71.043             304.630      -5.915     -35.139     -18.005     -33.177             0.000      -6.481       3.494              0.000            -13.948     0.000       7.351      24.518       0.000         0.230 1qys_0001\n",
      "\n",
      "2,3c2,3\n",
      "< SCORE: total_score       score dslf_fa13    fa_atr    fa_dun   fa_elec fa_intra_rep fa_intra_sol_xover4              fa_rep              fa_sol hbond_bb_sc hbond_lr_bb    hbond_sc hbond_sr_bb linear_chainbreak lk_ball_wtd       omega overlap_chainbreak            p_aa_pp pro_close rama_prepro         ref        time yhh_planarity description \n",
      "< SCORE:    -224.872    -224.872     0.000  -501.794   104.484  -146.069        1.118              18.786              71.043             304.630      -5.915     -35.139     -18.005     -33.177             0.000      -6.481       3.494              0.000            -13.948     0.000       7.351      24.518       0.000         0.230 1qys_0001\n",
      "---\n",
      "> SCORE: total_score       score dslf_fa13    fa_atr    fa_dun   fa_elec fa_intra_rep       fa_rep       fa_sol hbond_bb_sc hbond_lr_bb    hbond_sc hbond_sr_bb linear_chainbreak             omega overlap_chainbreak            p_aa_pp pro_close      rama       ref      time yhh_planarity description \n",
      "> SCORE:    -163.023    -163.023     0.000  -423.638   109.662   -46.146        1.040       49.117      241.309      -3.934     -26.998     -11.234     -25.491             0.000             4.211              0.000            -13.603     0.000    -4.905   -12.643     3.000         0.230 1qys_0001\n",
      "当前目录切换回到 /home/regen/git_develop/rosetta_notebook/00intro\n"
     ]
    }
   ],
   "source": [
    "with chdir(tutorial_path):\n",
    "    !$rosetta_home/source/bin/score_jd2.linuxgccrelease @flag > /dev/null\n",
    "    !cat output_files/score.sc\n",
    "    !echo\n",
    "    !diff output_files/score.sc output_files/expected_output/score.sc\n",
    "    !rm output_files/score.sc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 分数文件的分析\n",
    "\n",
    "注意官方提供的旧版ref2015打分结果与现在的差异.现在的分数是kcal/mol单位下的近似值.\n",
    "\n",
    "`total_score`表示总的加权分数.与未优化的结构相比,分数显著下降.\n",
    "\n",
    "对于这种大小的蛋白,通常分数在-300~-100REU.分数越低,给定蛋白的结构可能越稳定.对于ref2015,稳定结构每个残基提供约-3~-1REU.\n",
    "\n",
    "其他项是经过加权的结果,如`fa_atr`是加权后的残基间LJ吸引作用.这种分解可以帮助我们判断哪些能量项的贡献占主要部分.对于1QYS,残基间vdW作用对稳定的贡献最大,而溶剂作用对不稳定的贡献最大.\n",
    "\n",
    "`fa_rep`过大(例如比`fa_atr`的相反数大)说明结构中存在冲突.优化前后`fa_rep`明显下降了.\n",
    "\n",
    "### 更多打分选项\n",
    "\n",
    "#### 修改打分函数\n",
    "\n",
    "本节使用适合于多个蛋白对接docking的打分函数(具体权重可见`<path_to_Rosetta_directory>/main/database/scoring/weights`),使用的flag文件为flag_docking:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cd /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "当前目录从 /home/regen/git_develop/rosetta_notebook/00intro \n",
      "临时切换到 /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "-in:file:l input_files/pdblist\n",
      "\n",
      "-score:weights docking\n",
      "\n",
      "-out:file:scorefile output_files/score_docking.sc\n",
      "当前目录切换回到 /home/regen/git_develop/rosetta_notebook/00intro\n"
     ]
    }
   ],
   "source": [
    "with chdir(tutorial_path):\n",
    "    !cat flag_docking"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "执行打分程序:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cd /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "当前目录从 /home/regen/git_develop/rosetta_notebook/00intro \n",
      "临时切换到 /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "SEQUENCE: \n",
      "SCORE: total_score       score dslf_ca_dih dslf_cs_ang dslf_ss_dih dslf_ss_dst      fa_atr      fa_dun     fa_elec     fa_pair      fa_rep      fa_sol hbond_bb_sc hbond_lr_bb    hbond_sc hbond_sr_bb linear_chainbreak overlap_chainbreak               time description \n",
      "SCORE:    -111.103    -111.103       0.000       0.000       0.000       0.000    -169.606       5.373      -3.798      -2.237       5.683      73.720      -0.822      -8.609      -2.680      -8.128             0.000              0.000              0.000 1qys_0001\n",
      "SCORE:     -83.850     -83.850       0.000       0.000       0.000       0.000    -131.980       4.479      -3.295      -2.059       3.845      58.363      -1.663      -5.786      -1.092      -4.661             0.000              0.000              0.000 1ubq_0001\n",
      "当前目录切换回到 /home/regen/git_develop/rosetta_notebook/00intro\n"
     ]
    }
   ],
   "source": [
    "with chdir(tutorial_path):\n",
    "    !$rosetta_home/source/bin/score_jd2.linuxgccrelease @flag_docking > /dev/null\n",
    "    !cat output_files/score_docking.sc\n",
    "    # !rm output_files/score_docking.sc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "第1行为1qys的分数,第2行为1ubq的分数.此处的1qys总分数与孤立结构的分数不同.再次提醒,不同的打分函数的分数没有可比性.\n",
    "\n",
    "`dslf_fa13`能量项消失,变为4个描述双硫键几何性质的能量项.\n",
    "\n",
    "`fa_atr`分数有所变化,虽然1qys的结构以及vdW作用的计算公式是相同的,但权重不同.\n",
    "\n",
    "1ubq的总分数比1qys高,这是可预见的,因为rosetta是逐残基计算分数的,不能解释为1qys比1ubq稳定.目前不同蛋白的稳定性与总分数没有良好的相关性.\n",
    "\n",
    "#### 补丁文件,修改能量项权重\n",
    "\n",
    "假设我们想要修改docking打分函数中的权重,可能希望降低某些通用分数项的权重,提高其他分数项的权重,或者为约束添加分数项.有3种方式:\n",
    "\n",
    "- 提供自定义的权重文件,然后用`-score:weights`指定权重文件路径\n",
    "- 添加补丁文件,调整现有权重\n",
    "- 从命令行设置特定能量项的权重\n",
    "\n",
    "以第2/3种方式为例,我们提供以下`docking.wts_patch`补丁文件,对docking权重进行修正.例如将`fa_atr`从0.338修正为0.338*0.423=0.143；然后在命令行选项中将`fa_atr`权重重新设置为1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cd /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "当前目录从 /home/regen/git_develop/rosetta_notebook/00intro \n",
      "临时切换到 /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "fa_atr *= 0.423\n",
      "fa_rep *= 0.100\n",
      "fa_sol *= 0.372\n",
      "fa_intra_rep *= 0.000\n",
      "fa_pair *= 0.000\n",
      "fa_dun *= 0.064\n",
      "hbond_lr_bb *= 0.245\n",
      "hbond_sr_bb *= 0.245\n",
      "hbond_bb_sc *= 0.245\n",
      "hbond_sc *= 0.245\n",
      "p_aa_pp *= 0.00\n",
      "fa_elec = 0.026\n",
      "dslf_ss_dst *= 1.0\n",
      "dslf_cs_ang *= 1.0\n",
      "dslf_ss_dih *= 1.0\n",
      "dslf_ca_dih *= 1.0\n",
      "pro_close *= 0.000\n",
      "\n",
      "-in:file:l input_files/pdblist\n",
      "\n",
      "-score:weights docking\n",
      "-score:patch docking\n",
      "-score:set_weights fa_atr 1\n",
      "\n",
      "-out:file:scorefile output_files/score_docking_patch.sc\n",
      "当前目录切换回到 /home/regen/git_develop/rosetta_notebook/00intro\n"
     ]
    }
   ],
   "source": [
    "with chdir(tutorial_path) as (tmp,origin):\n",
    "    !cat $rosetta_home/database/scoring/weights/docking.wts_patch\n",
    "    !echo\n",
    "    !cat flag_docking_patch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cd /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "当前目录从 /home/regen/git_develop/rosetta_notebook/00intro \n",
      "临时切换到 /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "2,4c2,4\n",
      "< SCORE: total_score       score dslf_ca_dih dslf_cs_ang dslf_ss_dih dslf_ss_dst      fa_atr      fa_dun     fa_elec      fa_rep      fa_sol hbond_bb_sc hbond_lr_bb    hbond_sc hbond_sr_bb linear_chainbreak overlap_chainbreak               time description \n",
      "< SCORE:    -482.214    -482.214       0.000       0.000       0.000       0.000    -501.794       0.344      -3.798       0.568      27.424      -0.201      -2.109      -0.657      -1.991             0.000              0.000              0.000 1qys_0001\n",
      "< SCORE:    -374.621    -374.621       0.000       0.000       0.000       0.000    -390.474       0.287      -3.295       0.384      21.711      -0.407      -1.418      -0.268      -1.142             0.000              0.000              0.000 1ubq_0001\n",
      "---\n",
      "> SCORE: total_score       score dslf_ca_dih dslf_cs_ang dslf_ss_dih dslf_ss_dst      fa_atr      fa_dun     fa_elec     fa_pair      fa_rep      fa_sol hbond_bb_sc hbond_lr_bb    hbond_sc hbond_sr_bb linear_chainbreak overlap_chainbreak               time description \n",
      "> SCORE:    -111.103    -111.103       0.000       0.000       0.000       0.000    -169.606       5.373      -3.798      -2.237       5.683      73.720      -0.822      -8.609      -2.680      -8.128             0.000              0.000              0.000 1qys_0001\n",
      "> SCORE:     -83.850     -83.850       0.000       0.000       0.000       0.000    -131.980       4.479      -3.295      -2.059       3.845      58.363      -1.663      -5.786      -1.092      -4.661             0.000              0.000              0.000 1ubq_0001\n",
      "当前目录切换回到 /home/regen/git_develop/rosetta_notebook/00intro\n"
     ]
    }
   ],
   "source": [
    "with chdir(tutorial_path):\n",
    "    !$rosetta_home/source/bin/score_jd2.linuxgccrelease @flag_docking_patch > /dev/null\n",
    "    !diff output_files/score_docking_patch.sc output_files/score_docking.sc\n",
    "    !rm output_files/score_docking*.sc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看到不同项之间的比例与补丁文件和命令行设置对应.注意由于将`fa_pair`权重比例设置为0,新的打分文件中不会出现该列.\n",
    "\n",
    "#### 高级选项\n",
    "\n",
    "其他flag文件相关的选项参见[打分命令](#TODO).\n",
    "\n",
    "例如`-scorefile_format json`将输出文件格式设置为json.\n",
    "\n",
    "### 获得单个残基的分数\n",
    "\n",
    "`total_score`是单个残基分数的总和.获取每个残基的分数对于详细分析很有帮助,例如对于找出有空间冲突的残基特别有用.\n",
    "\n",
    "本节使用`per_residue_energies`程序,输出格式为silent(`-out:file:silent`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cd /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "当前目录从 /home/regen/git_develop/rosetta_notebook/00intro \n",
      "临时切换到 /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "-in:file:s input_files/1qys.pdb\n",
      "\n",
      "-out:file:silent output_files/per_res.sc\n",
      "\n",
      "SCORE:  pose_id               pdb_id  fa_atr  fa_rep  fa_sol  fa_intra_rep  fa_intra_sol_xover4  lk_ball_wtd  fa_elec  pro_close  hbond_sr_bb  hbond_lr_bb  hbond_bb_sc  hbond_sc  dslf_fa13  omega   fa_dun  p_aa_pp  yhh_planarity  ref     rama_prepro  score   description\n",
      "SCORE:  input_files/1qys.pdb  3A      -2.902  0.187   3.047   0.011         0.940                -0.046       -1.618   0.000      0.000        0.000        0.000        -0.828    0.000      0.019   1.594   0.000    0.000          -2.146  0.000        -1.743  residue_1\n",
      "SCORE:  input_files/1qys.pdb  4A      -7.106  0.617   2.317   0.034         0.096                -0.041       -1.661   0.000      0.000        -0.734       0.000        0.000     0.000      0.001   0.785   -0.448   0.000          2.304   -0.160       -3.997  residue_2\n",
      "SCORE:  input_files/1qys.pdb  5A      -5.601  1.002   5.613   0.009         0.500                0.427        -3.629   0.000      0.000        -1.490       0.000        -1.216    0.000      -0.002  3.033   0.281    0.000          -1.451  -0.068       -2.592  residue_3\n",
      "SCORE:  input_files/1qys.pdb  6A      -6.259  1.078   1.135   0.018         0.048                0.081        -1.953   0.000      0.000        -1.295       0.000        0.000     0.000      0.033   0.032   -0.594   0.000          2.643   -0.047       -5.081  residue_4\n",
      "当前目录切换回到 /home/regen/git_develop/rosetta_notebook/00intro\n"
     ]
    }
   ],
   "source": [
    "with chdir(tutorial_path):\n",
    "    !cat flag_per_residue\n",
    "    !echo\n",
    "    !$rosetta_home/source/bin/per_residue_energies.linuxgccrelease @flag_per_residue > /dev/null\n",
    "    !head -n 5 output_files/per_res.sc | column -t\n",
    "    !rm output_files/per_res.sc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "pdbid类似`3A`,表示A链的第3个残基,3是pdb编号；而最后一列描述为residue_i,i从1开始,是rosetta编号.\n",
    "\n",
    "输出结果并没有对齐,因此使用column命令进行对齐.也可以使用Excel,Calc等办公软件打开查看.如果要将`per_residue_energies`和`score_jd2`的结果进行对比,确保使用相同的打分函数.\n",
    "\n",
    "### 获得单个残基的分数明细\n",
    "\n",
    "如果要进一步分解残基的能量,可以使用`residue_energy_breakdown`程序获得残基自身和其他残基的相互作用分数.flag文件和打分结果如下:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cd /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "当前目录从 /home/regen/git_develop/rosetta_notebook/00intro \n",
      "临时切换到 /home/regen/git_develop/rosetta/demos/tutorials/scoring\n",
      "-in:file:s input_files/1qys.pdb\n",
      "\n",
      "-out:file:silent output_files/energy_breakdown.sc\n",
      "\n",
      "SCORE:  pose_id               resi1  pdbid1  restype1  resi2  pdbid2  restype2  fa_atr  fa_rep  fa_sol  fa_intra_rep  fa_intra_sol_xover4  lk_ball_wtd  fa_elec  pro_close  hbond_sr_bb  hbond_lr_bb  hbond_bb_sc  hbond_sc  dslf_fa13  omega   fa_dun  p_aa_pp  yhh_planarity  ref     rama_prepro  total   description\n",
      "SCORE:  input_files/1qys.pdb  1      3A      ASP       --     --      onebody   0.000   0.000   0.000   0.011         0.940                0.000        0.000    0.000      0.000        0.000        0.000        0.000     0.000      0.019   1.594   0.000    0.000          -2.146  0.000        0.418   input_files/1qys.pdb_1_onebody\n",
      "SCORE:  input_files/1qys.pdb  1      3A      ASP       2      4A      ILE       -1.671  0.105   1.404   0.000         0.000                -0.199       0.076    0.000      0.000        0.000        0.000        0.000     0.000      0.000   0.000   0.000    0.000          0.000   0.000        -0.285  input_files/1qys.pdb_1_2\n",
      "SCORE:  input_files/1qys.pdb  1      3A      ASP       3      5A      GLN       -1.820  0.087   3.264   0.000         0.000                0.627        -3.146   0.000      0.000        0.000        0.000        -1.657    0.000      0.000   0.000   0.000    0.000          0.000   0.000        -2.645  input_files/1qys.pdb_1_3\n",
      "SCORE:  input_files/1qys.pdb  1      3A      ASP       20     22A     THR       -0.957  0.164   0.685   0.000         0.000                -0.313       -0.082   0.000      0.000        0.000        0.000        0.000     0.000      0.000   0.000   0.000    0.000          0.000   0.000        -0.504  input_files/1qys.pdb_1_20\n",
      "当前目录切换回到 /home/regen/git_develop/rosetta_notebook/00intro\n"
     ]
    }
   ],
   "source": [
    "with chdir(tutorial_path):\n",
    "    !cat flag_residue_energy_breakdown\n",
    "    !echo\n",
    "    !$rosetta_home/source/bin/residue_energy_breakdown.linuxgccrelease @flag_residue_energy_breakdown > /dev/null\n",
    "    !column -t output_files/energy_breakdown.sc |grep -E 'pose_id| 3A ' |head -n 5\n",
    "    !rm output_files/energy_breakdown.sc*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "标题行后的第1行是3A残基与自身的相互作用,后续有3A和其他残基的相互作用分数.\n",
    "\n",
    "### 高级打分流程\n",
    "\n",
    "有些蛋白需要额外的打分文件.例如[膜蛋白](#TODO)需要额外的膜区域定义文件,[对称蛋白](#TODO)需要对称性的定义文件.具体流程见对应链接.\n",
    "\n",
    "## 参考内容\n",
    "\n",
    "- 对rosetta打分函数的总结文章:[Alford et al. (2017) J. Chem. Theory Comput. In press.](https://www.ncbi.nlm.nih.gov/pubmed/28430426)\n",
    "- ref2015实现与性能文章:[Park et al. (2016). J. Chem. Theory Comput. 12(12):6201-12](https://www.ncbi.nlm.nih.gov/pubmed/27766851)\n",
    "- talaris2014实现与性能:[O’Meara et al., J. Chem. Theory Comput. 2015](https://dx.doi.org/10.1021/ct500864r)\n",
    "- [talaris2013的修改内容](https://www.rosettacommons.org/node/3508#comment-6946)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
